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Abstract 

In the mean-field approximation, the number density p(r) for the ground state of 
a Bose-Einstein condensate trapped by an external potential V{r) satisfies a classical 
field equation called the Gross-Pitaevskii equation. We show that quantum corrections 
to p are dominated by quantum fluctuations with wavelengths of order where 
a is the S-wave scattering length. By expanding the equations for the Hartree-Fock 
approximation to second order in the gradient expansion, we derive local correction 
terms to the Gross-Pitaevskii equation that take into account the dominant effects of 
quantum fluctuations. We also show that the gradient expansion for the density breaks 
down at fourth order. 
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1 Introduction 



The successful achievement of Bose-Einstein condensation of atomic gases in magnetic traps 
|1|, 0, ^ has created an explosion of interest in interacting Bose gases. The condensates 
in existing magnetic traps are sufficiently dilute that the mean- field approximation gives a 
satisfactory description of present experimental measurements. However accurate theoretical 
predictions require that quantum fiuctuations around the mean field also be taken into 
account. The relative magnitude of these corrections grows as the square root of the number 
density of the atoms. They will therefore become more important as higher trap densities 
are achieved and as the precision of experimental measurements improves. 

One of the basic observables of a Bose-Einstein condensate trapped in an external poten- 
tial V{r) is the number-density profile p(r) of the ground state. In the mean-field approxi- 
mation, p(r) satisfies the Gross-Pitaevskii equation: 

+ - ^(r) j v^(r) - ^^pM^) = ' (1) 

where a is the S-wave scattering length of the atoms. The chemical potential must be 
tuned so that / d^r p = N, where N is the number of atoms in the trap. The density profile 
of a trapped Bose-Einstein condensate has been studied extensively using (p. The solutions 
to this equation have been calculated using numerical methods |^ and variational methods 
B 0' H' El- '^^^ solutions have also been studied analytically in the Thomas- Fermi limit, in 
which the gradient term in (|l|) is neglected [§] . The corrections due to the breakdown of this 



approximation near the edge of the condensate have also been studied ^ . 

There are corrections to the mean-field approximation from quantum fiuctuations around 
the mean field. In a dilute homogeneous Bose gas, the relative magnitude of the contributions 
of quantum fiuctuations to thermodynamic quantities is characterized by the dimensionless 
quantity \/pa^. For condensates in existing magnetic traps, the peak value of \/pa^ is small 
but not negligible. Since there are some observables that are more sensitive than the density 
to the effects of quantum fluctuations, it is important to be able to calculate the effects of 
quantum fluctuations quantitatively. 

In this paper, we calculate the effects of quantum fluctuations on the density proflle for 
a Bose-Einstein condensate in a trapping potential. The expansion parameter that charac- 
terizes the relative magnitude of these effects is y/pa^-, where p is the local number density. 
We point out that the quantum corrections are dominated by quantum fluctuations with 
wavelengths on the order of 1/^Jpa. The leading effects of these short-distance quantum 
fluctuations can be calculated using the gradient expansion. By carrying out a self-consistent 
one-loop calculation through second order in the gradient expansion, we determine the cor- 
rection terms that must be added to the Gross-Pitaevskii equation to take into account 
the effects of quantum fluctuations: 
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Our method involves a combination of the Hartree-Fock approach [12| and the Thomas- Fermi 
approach |jl3|, In the Hartree-Fock approximation, which involves the self-consistent 
treatment of one-loop quantum corrections, the equation for the density is an integral equa- 
tion. We obtain the local differential equation by applying a gradient expansion to the 
integral equation, which corresponds to expanding around the Thomas- Fermi limit. 

We begin in Section ^ by formulating the problem of Bose-Einstein condensation in a 
trapping potential as a problem in quantum field theory. In Section ^ we develop a per- 
turbation expansion for calculating the effects of quantum fluctuations around an arbitrary 
background field. In Section ^ we calculate the one-loop quantum corrections to the density 
profile and the condensate profile. We show that the ultraviolet divergences that arise in the 
calculation can be removed by the same renormalizations of the action and the number den- 
sity that are required in the absence of the potential. We find that the number density can 
be expanded to second order in the gradient of the mean field, while the gradient expansion 
for the condensate breaks down at that order. In Section ^, we calculate the self-consistent 
one-loop quantum corrections to the density profile to second order in the gradient expan- 
sion and show that they are given by (^. We repeat the calculation in Section |^ using an 
alternative parameterization for the quantum field and demonstrate that the final result is 
independent of the parameterization. We also use this parameterization to show that the 
gradient expansion for the density breaks down at fourth order. Finally, in Section ^ we 
examine the implications of our results for Bose-Einstein condensation in present magnetic 
traps. Details of the calculations of Feynman diagrams are included in several appendices. 



2 Quantum field theory formulation 



Consider a large number N of identical bosonic atoms trapped in an external potential V{r) 
at zero temperature. If the momenta of the atoms are sufficiently low, their de Broglie wave- 
lengths are much smaller than the range of the interatomic potential, which is comparable 
in magnitude to the Bohr radius. In this case, the scattering of two atoms will be dominated 
by S-wave scattering and can be described by a single number, the S-wave scattering length 
a. Our problem is to determine the number- density profile p(r) of the atoms. We begin 
by formulating this many-body quantum mechanics problem as a problem in quantum field 
theory. 

A convenient way to describe a system containing any number of atoms is in terms of 
a quantum field ip{r,t) that annihilates an atom. If the atom is a boson, the field satisfies 
the equal-time commutation relations 



V(r,t),V^(r',t) 
V'(r,t),^t(r',t) 



0, 
5{y 



The time evolution of the field is given by the equation 



(3) 
(4) 



(5) 
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where the couphng constant g is related to the S-wave scattering length a by 



9 



m 



(6) 



The parameter Sg in is a counterterm associated with the renormalization of g. It is 
necessary to impose an ultraviolet cutoff Auv on the wavenumbers of virtual particles in order 
to avoid ultraviolet divergences due to quantum fluctuations of the field. Renormalization 
of a quantum field theory is the tuning of its parameters so that physical quantities are 
independent of the ultraviolet cutoff. All the dependence of first order quantum corrections 
on Auv can be removed by adjusting 6g in as a function of Auy. 
The number operator, which counts the number of atoms, is 

U = j d^rij^^lj{r,t) . (7) 

That this is a number operator follows from the commutation relations @ and @, which 
imply that ip'^ and if) act as raising and lowering operators for M. The equation implies 
that Af is independent of time, so the number of atoms is conserved. The equation can 
also be expressed in the form 

ih^i^ = -[n,^p] , (8) 



dt 



where the hamiltonian operator Ti is 



n 



9 + ^9 



(9) 



The hamiltonian Ti measures the energy of the atoms. 

The vacuum |0), defined by '?/'(r, t)|0) = 0, represents the state containing zero atoms. 
One can show that a Schroedinger wavefunction for N atoms can be represented as a matrix 
element of an operator between a state with Af = N and the vacuum. The simplest case is 
a state |0) containing 1 atom, which satisfies A/'|0) = \(f))- Since the last term in @ annihi- 
lates the single particle state the matrix element {0\tp{r, t)\(f)) satisfies the Schroedinger 
equation: 

d 



ih 



dt 



Viv) 



(O|V^(r,t)|0) 



0. 



(10) 



Thus {0\ip{r,t)\(f)) can be interpreted as the Schroedinger wavefunction for an atom in the 
potential ^(r). 

The next simplest case is a state |0) containing 2 atoms, which satisfies Af\(j)) = 2|0). It is 
straightforward to show using that the matrix element {0\ip{ri,t)ip{r2,t)\(j)) satisfies the 
Schroedinger equation for two particles in the external potential V"(r) interacting through 
a two-body potential proportional to 6^{ri — r2). In the absence of the potential V^(r), one 
can calculate the amplitude for the scattering of two atoms exactly |jl5|. The scattering 
amplitude f{9) is independent of the scattering angle 9, so it describes S-wave scattering. If 
the total energy of the two atoms in the center-of-momentum frame is £" = 2(p^/2m), the 
scattering amplitude is 



/ 



1 

47r 



2h' 



m{g + 6g) 
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d^k 



{2n)^ k^ -mE/n' -ie 



-1 



(11) 
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The integral over the wave- vector k is ultraviolet divergent. A particularly convenient reg- 
ularization of the integral is to introduce an ultraviolet cutoff Auv through the following 
prescription: 

d^k 1 1 /"^ ,9 „ / 1 1 \ 



(27r)3 P -m^/r -ie 2n^ Jo \k^ -mE/fi^ ~ie P 



1 f^uv „ 1 



The scattering amplitude then becomes 



^ Air 



2%^ 1 . 

+ —Auv + 



m[g + 5g) 271^ A'nh 



(13) 



The dependence on the ultraviolet cutoff can be completely cancelled by choosing the bare 
coupling constant g + 5g to be 

^ + ^^ = 1 - mgL/{2r:hf ^''^ 
The scattering amplitude (|T3|) then reduces to 



This confirms the identification of a as the S-wave scattering length. The scattering of atoms 
is correctly reproduced by the pointlike interaction term in as long as the energy of the 
atoms is sufficiently low that (0) is a good approximation to the scattering amplitude. Note 
that the energy dependence in (p!5D is that required by the optical theorem. 

It is sometimes stated that a delta-function potential in 3 dimensions is trivial in the 
sense that it gives no scattering. A more accurate statement is that there is no scattering if 
we take the ultraviolet cutoff to infinity with the strength of the potential held fixed. This is 
evident from (|I3|), which shows that / ^ if we take Auv oo with g + Sg fixed. However, 
if we allow the strength of the potential to vary with Auv in accordance with (p!4D, we obtain 
the nontrivial scattering amplitude (]T3|). 

From the expression (p!3D for the scattering amplitude, one can infer an upper limit on 
the ultraviolet cutoff that must be satisfied in order for perturbative calculations to be 
accurate. The expansion for / in powers of g, including the first quantum correction which 
is proportional to g"^, is 



^ mg 



( mgAijy mgVmE\ 



(16) 



If the ultraviolet cutoff Auv is too large, there is a delicate cancellation between the term 
proportional to Auv in (P^, which comes from an integral over k, and the counterterm 6g. 
Since a perturbative calculation is necessarily approximate, the cancellation can lead to large 
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errors. Such a large cancellation can be avoided if the term proportional to Auv in (0) is 
much less than 1. This sets an upper bound on the ultraviolet cutoff: 



Auv « ^ = ^. 17 

mg 2a 

If this upper bound on Auv is not satisfied, then in order to obtain an accurate calculation, 
it is necessary to use a nonperturbative calculational method that sums up all orders in g. 

A state 10) containing 3 atoms satisfies A^l^) = 3|0). In the absence of the potential 
V{r), one can calculate the amplitude for 3 — > 3 scattering as an expansion in powers of 
g. The leading contribution is proportional to g^ and comes from two successive 2 —>■ 2 
scatterings. Higher order terms in g represent quantum corrections. The dependence of the 
first quantum correction on the ultraviolet cutoff Auv is cancelled by the counterterm 6g in 
(^). However, the second quantum correction, which is proportional to g'^ has a logarithmic 
ultraviolet divergence that is not cancelled [|r^ . Thus corrections to physical quantities from 
second order in the quantum fluctuations depend on the ultraviolet cutoff. One can eliminate 
the dependence on Auv from second order quantum corrections by adding to (||) the term 
((73 + Sg3)'il:^ip^ipil'xl'/18. The logarithmic ultraviolet divergence is cancelled by choosing the 
counterterm to be 

S93 = ^TTl ^ 9 log , (18) 

where the value of n, which was introduced to make the argument of the logarithm dimen- 
sionless, depends on the precise definition of g^. The parameter g^ represents a point-like 
contribution to the 3^3 scattering amplitude. Thus the S-wave scattering length a does 
not contain enough information about the low-energy scattering of atoms to calculate second 
order quantum fluctuations. It is also necessary to know the 3 — > 3 coupling constant g^. In 
this paper, we will avoid this complication by calculating only to first order in the quantum 
corrections. 

A state 10) containing atoms satisfies A/'|0) = A^|0). In the presence of the potential 
V{r), the ground state in the Af = N sector, which we denote by I^at), is the state that 
minimizes (0|7i|0) subject to the constraint Af\(p) = N\(f)). The desired number-density 
profile is 

p(r) = (fi^lV^V(r)lfi^). (19) 

If N is large, we expect p(r) to be insensitive to changes in N that are small compared 
to N. This suggests that we can relax the constraint on the particle number and replace 
the state \ flN) in (IHD by the state that minimizes (0|7i|0) subject to the weaker constraint 
(0|A/'|0) = A^. If the root-mean-square fluctuations of A/" in that state are small compared 
to A^, the expectation value of ip'^ip^r) in that state should give a good approximation to 
(|I^). But that state is precisely the ground state \ fl^) of the quantum field theory whose 
hamiltonian is ?i — fiAf, where the value of the chemical potential n is such that 

{n^\Af\n,) = N. (20) 

Thus, if A^ is sufficiently large, the density profile can be approximated by the ground-state 
expectation value of the operator ip'^ip^r) in the state \^^)- 
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We have now formulated the problem of calculating the density profile as a quantum field 
theory problem. The field theory is summarized by the action 



S[^/j] = dt d^rl 



^-^(^Vn. (21) 



The counterterms 5/i and 5g are needed to cancel ultraviolet divergences associated with 
quantum fluctuations of the field. The counterterm 5/1 would also have been required in (^) 
if the ip^ipip term had not been normal-ordered. The local number- density operator ■?/'"'"■?/; (r) 
can also have ultraviolet divergences, and a counterterm 6p is therefore needed to cancel 
these divergences. The number-density profile is the expectation value of this local operator 
in the ground state of the field theory: 

p(r) = (^V(r)) - 6p. (22) 

Here and below, we use the angular brackets (...) to denote the expectation value in the 
ground state The chemical potential /z in ( pT]) must be adjusted so that the integral 

of the local number density is equal to the number of atoms: 

J d\ p{r) = N . (23) 

Thus our problem reduces to calculating the ground-state expectation value (E2[) for the 



quantum field theory defined by (|2l|) . Another important observable is the condensate 
{ip{r)), which is the ground-state expectation value of the field. A nonzero value of the 
condensate indicates the spontaneous breaking of the phase symmetry ip e^'^ip of the 
action (PT]). 

The ultraviolet divergences that are cancelled by the counterterms 5g, and 5p arise 
from treating the interaction between atoms as pointlike down to arbitrarily short distances. 
The divergences could be avoided by replacing the pointlike {ip'^ipY interaction term in the 
action (|2l|) by an interaction through a two-body potential v{jc 1 — ^2)- A physically realistic 
two-body potential would have a range comparable to the size of an atom and its shape 
would have to be adjusted so that it gives the correct S-wave scattering length a. It would 
be rather inefficient to calculate the effects of interactions using a physically realistic two- 
body potential. The reason is that physical quantities depend on the two-body potential 
in a very simple way. Almost all of the dependence enters through the S-wave scattering 
length a. Thus we can obtain the same result for physical quantities by using a pointlike 
interaction that gives the same scattering length. A pointlike interaction necessarily leads to 
ultraviolet divergences from short-distance quantum fluctuations, and they must be cancelled 
by counterterms. Since the counterterms cancel short-distance effects, they correspond to 
renormalization of the parameters in the action and to renormalization of local operators, 
such as ip'^ipir). We will find that the three counterterms 6g, and 6p are sufficient 
to cancel all the ultraviolet divergences associated with first order quantum fluctuations. 
The values of the counterterms are independent of the external potential V{r). Thus all 
ultraviolet divergences are removed by the same renormalizations that are required for the 
homogeneous Bose gas. 

An alternative way to deal with ultraviolet divergences is to replace the interatomic 
potential {g + 6g)S^{r)/2 by a pseudopotential g6^{r){d/dr)r/2 Ultraviolet divergences 



7 



can be avoided by evaluating the partial derivative in the pseudopotential at the appropriate 
stage of the calculation. We find it simpler to introduce an ultraviolet cutoff and use the 
renormalization machinery of quantum field theory to remove the dependence on the cutoff. 

Most previous work on the density profile has been carried out within the mean-field 
approximation. In this approximation, quantum fiuctuations are neglected. The countert- 
erms 6fi, 6g, and 6p, which cancel ultraviolet divergences associated with those quantum 
fiuctuations, can all be set to zero. The field ijj satisfies the time-dependent Gross-Pitaevskii 
equation, which is the classical field equation for the action (|2T|): 

'^Ft"^ + (1^^' + -"-^^)^ - f (^V)^ = 0. (24) 

The ground state \ fl^) corresponds to a time-independent solution 0o(r) that can be chosen 
to be real-valued. The mean field therefore satisfies 

(l^V^ + ^ - ^(r)) Mr) - f 0g(r) = . (25) 
The number density (^) reduces to 

p(r) = 0^(r) . (26) 

Thus the density profile in the mean-field approximation satisfies the Gross-Pitaevskii equa- 
tion (|l]). The condensate in the mean-field approximation is 

(^(r)) = 0o(r) , (27) 

and it therefore satisfies (ip) = y/p. 

The density profile is modified by quantum corrections. The corrections can be obtained 
by expanding the quantum field around the mean field (po{r), which satisfies ([25|): 

ijir,t) = 0o(r) + ^(r,t). (28) 

The expression (^2|) for the number density becomes 

p(r) = (Plir) + 20o(r)Re(V^(r)) + (^t^(r)) - 6p. (29) 

The condensate differs from y/p because of the effects of quantum fluctuations: 

(V'(r)) = Mr) + (^(r)) . (30) 

Having formulated our problem in terms of quantum fleld theory, there are quantum 
fluctuations on all length scales ranging from L, the length scale associated with variations 
in p(r), down to the inverse of the ultraviolet cutoff Ayy. The quantum fluctuations with 
length scales of order L depend in detail on the shape of the potential V^(r). For quantum 
fluctuations with length scales much smaller than L, the effects of variations in V{r) are 
negligible. The short-distance quantum fluctuations therefore behave locally like those of a 
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homogeneous Bose gas with chemical potential fi — V{r). We will show that these fluctuations 
give the dominant corrections to the density proflle. 

The one-loop quantum corrections to the number density can be obtained by keeping 
the terms in the action that are quadratic in the fluctuating flelds ilj{r,t). If these flelds 
are expanded in terms of normal modes, the corrections ( pOf ) to p(r) can be expressed as a 
sum over the normal modes. The contribution of an individual normal mode to the number 
density scales like This is negligible compared to the contribution from the mean 

fleld, which scales like N/L^. A signiflcant contribution can only arise from summing over 
a large number of normal modes. Normal modes with very short wavelengths approach a 
continuum and can be labeled by the wavevector k. The contribution to the density from 
such modes scales like / d^k. The integral is ultraviolet divergent. The ultraviolet divergence 
is proportional to A^y and is cancelled by the counterterm 6p in p9|). After renormalization, 
modes with k comparable to Auv do not contribute to p(r). Since the density of modes grows 
rapidly with k, the dominant quantum corrections to p(r) come from the largest values of 
k whose effects are not removed by renormalization. To understand the scale of k that 
dominates, it is useful to recall some simple facts about the homogeneous Bose gas. 

The properties of a homogeneous Bose gas with positive scattering length a and low 
number density p are well understood. The dimensionless quantity i/po^ serves as an expan- 
sion parameter for the low-density expansion. For example, the ground-state energy density, 
including the flrst quantum correction, is 



S = —^1 + TT^Vpa3 . (31) 
The coefficient of \/pa^ in the quantum correction term was flrst obtained by Lee and Yang 



18| . The quasiparticle excitations of the system are Bogoliubov modes, which are plane 



waves with the dispersion relation 

. 5!^-^. (32) 
This dispersion relation changes from linear in k to quadratic at a scale A given by 



A = ^IQTrap. (33) 

This is the scale of the wavenumber k that dominates the quantum corrections to the energy 
density. The one-loop quantum correction is the sum over normal modes of the zero-point 
energies huj/2, where u is the angular frequency of the normal mode. The contribution 
from large k behaves like / d^ke{k)/2, where e{k) is the Bogoliubov dispersion relation given 
in (|32D. This integral is ultraviolet divergent, with the leading divergence proportional to 
fi^A^y/m. This leading divergence and the subleading divergences can all be removed by 
renormalization. After renormalization, the integral is dominated by the scale A given in 
(^), and it therefore scales like h^A^/m /m. This estimate agrees with the 

explicit result given in (pl|). 

Generalizing to the case of a nonhomogeneous Bose gas, we can anticipate that the 
quantum corrections to the density proflle p(r) will be dominated locally by modes with 
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wavenumber k on the order of Yl6'7rap(r). The contributions from much shorter wave- 
lengths are removed by renormahzation. The contributions from much longer wavelengths 
are suppressed by phase space. These modes can be approximated by a continuum as long 
as the corresponding wavelengths are much shorter than the scale L for significant variations 
in p(r): 

y^l67rap(r) > ^. (34) 

If this lower bound on the density is satisfied, then the methods of continuum quantum field 
theory can be used to calculate the dominant quantum corrections. The condition (|3^ ) is 
also necessary in order to calculate quantum corrections using a gradient expansion, which 
is an expansion in 1/ (AL). There is an upper bound on p(r) that must be satisfied in order 
to allow perturbative calculations in the quantum field theory with the pointlike interaction 
in (pri). This condition is that the scale A must be much less than the maximum ultraviolet 
cutoff given in (p!7|): 

\/l6vrap(r) « ^. (35) 

If this upper bound is not satisfied, then nonperturbative methods must be used to calculate 
the effects of quantum corrections. The condition (^) is also necessary in order for the 
effects of quantum fluctuations to be small enough to be treated as perturbative corrections 
to the mean-field approximation. For example, in the case of a homogeneous Bose gas, the 
condition that the quantum correction to the energy density, which is given by the second 
term in (PTD, is small compared to the mean-field contribution is essentially identical to (|53|). 
In our analysis of quantum corrections to the density profile, we will assume that the number 
density is in the range specified by the inequalities (0) and (pSl) . 



3 Perturbative Framework 

In this Section, we present a general framework for carrying out perturbative calculations of 
the effects of quantum fluctuations around an arbitrary time-independent background v{y). 
In Section ^ we will set v equal to the mean field 0o and use this framework to calculate 
one-loop corrections to the condensate and density profiles. In Section ^, we set v equal to 
the condensate {ip) and determine the self-consistent one-loop corrections to the equation 
for the density profile. 

It is convenient to parameterize the quantum field t) in terms of two real- valued quan- 
tum fields ^ and rj that describe quantum fluctuations around an arbitrary time-independent 
background v{r): 

iP{r,t) = v{r) + . (36) 

We will refer to this as the Cartesian parameterization of the quantum field. An alternative 
field parameterization is considered in Section If the phase of ip is chosen so that v is 
real-valued, the condensate profile is 

(^(r)) = v{r) + ^(eW). (37) 
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The number density is 
p(r) = V 



\r) + V2v{r)m) + \ {e{r)) + \ {ri\v)) - Sp. 



(38) 



Inserting the field parameterization (|36D into the action (^) and expanding in powers of 
the quantum fields ^ and r], it becomes 



S[v] + JdtJ f/vjv^Te +l{v^~ ^v) + (V^ - A' + X) e 



9 + 



16 



,(39) 



where / = |^/ and T, X, Y, and Z are external sources that depend on v: 



T(r) 
X(r) 
y(r) 



9 + ^9 2/ X 
(r) 



^(r) + — V^f(r) 
2m 



+ 2m 



2m 
9 + ^9 



-vir 



(40) 
(41) 
(42) 
(43) 



We have set ^ = 1 in the action. Dimensional analysis can be used to reinsert the factors of h 
at the end of the calculation. The parameter A appears both in the source X and explicitly 
in the term in the action, and cancels between them. The arbitrariness of this parameter 
can be exploited to simplify calculations. 

To organize the quantum corrections into a loop expansion, we separate the terms in the 
action that depend on ^ and rj into a free part and an interaction part: 



Slip] = S[v] + Sfree[^,r]] + Sint[v,^,r]] 

The free part of the action is 



dt / £r{ - (r]i - ^f]) + — ^(V^ - A^)^ + —r]V^r] 



Am 



Am 



(44) 



(45) 



This action describes Bogoliubov modes with the dispersion relation (p2D, where A is now 
an adjustable parameter. The Fourier transform of the propagator for the fields ^ and 7] is 
a 2 X 2 matrix: 



/ D««(k,cu) Z^«''(k,cu) 



1 



e/{2m) 



-luo 



(46) 



where k is the wavevector and uj is the frequency. The diagonal elements of the propagator 
matrix (|46|) are represented by solid lines for ^ and dashed lines for ry, as illustrated in 
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Figures 1(a) and 1(b). The off-diagonal elements are represented by a line that is half solid 
and half dotted, as in Figure 1(c). All the remaining terms in the action (|39[) are treated as 



interactions: 

Sintlf,?,-)] = /*/ d'r^^/2T( + + ^Y'f 

They include interactions with the external sources T, X, Y, and Z as well as the 4-point 
coupling g + 6g. The sources are represented by dots labeled by the appropriate letter, as 
illustrated in Figure 2. The 4-point couplings are represented by points that connect four 
lines. 

It is possible to diagonalize the propagator matrix (^6|) by applying a Bogoliubov trans- 
formation to the fields ^ and rj. However, such a transformation makes the interaction terms 
in the action significantly more complicated and increases the number of diagrams that 
contribute to most quantities. For explicit calculations, it is more economical to minimize 
the number of diagrams. We therefore prefer to use a propagator matrix with off-diagonal 
elements. 

In the case of a homogeneous Bose gas, the trapping potential V is zero and we can choose 
the background field f to be a constant, independent of r. If we choose f ^ = 2(yU+5/i) / (g+Sg) 
and = 2m{g + ^g)v'^, then the sources T, X, and Y in (^) vanish and the interactions 
reduce to three-point couplings and four-point couplings. Such a perturbative framework 
has recently been used to reproduce the classic one-loop corrections to the thermodynamic 



properties of a homogeneous Bose gas |19 



The leading quantum corrections to the ground-state expectation values in (pTf ) and ( p8D 
are given by one-loop Feynman diagrams. Examples of one-loop diagrams that contribute 
to {^^) and {t]^) are shown in Figures 3 and 4, respectively. The dot on the left side of 
each diagram represents the operator which creates two solid lines, or the operator rj'^, 
which creates two dashed lines. The lines form a loop that can include any number of 
insertions of the sources X and Y. It is convenient to introduce the notation (. . .)xy for the 
expectation value of an operator in the presence of the sources X and Y, but with no other 
self-interactions for the quantum fields. The sum of all one-loop diagrams for (^^) and {rj'^) 
can then be represented as 



e{r) 



r/2(r) 



1— loop 
1— loop 



('M),,. , (48) 



')'(r))^j. . (49) 



The advantage of this notation is that the expectation value {^) at one-loop order can also be 
expressed succinctly in terms of (^^)x,y and {rf)x,Y- Examples of diagrams that contribute 
to (^) are shown in Figures 5, 6, and 7. The operator ^ creates a single solid line. In the 
diagrams of Figure 5, the ^ field propagates to a source T, where it is annihilated. In the 
diagrams of Figures 6 and 7, it propagates to a source Z, which creates a pair of solid lines or 
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dashed lines that form a loop. In all of these diagrams, the ^ propagator and the propagators 
inside the loop can include any number of insertions of the sources X and Y. The sum of 
all such diagrams can be expressed as 



(e(r)) 



1— loop 



-v^ J d^r' (^j dt' D^l 



X (r(rO + \z{r'){e{v'))^^^ + \z{j'){ri\v'))^^^] , (50) 



where -D^ y is the diagonal component of the propagator for ^ in the presence of the sources 
X and Y.' 

The quantities (^^)x,y and {rf')x,Y in (^HD, (^9]), and ( pOD are functional of the sources 
X and Y . These functionals include terms of arbitrarily high orders in X and Y . They are 
nonlocal because the loop diagrams involve an integral over the positions of the sources X 
and Y . After renormalization, these integrals are dominated by wavelengths of order 27r/A, 
while the sources vary significantly only over much larger distances of order L. It is therefore 
reasonable to expand the sources X{y') and Y{y') as Taylor series around the point r. This 
reduces the expressions for (^^)x,y and {rf)x,Y to an infinite sum of local quantities involving 
X, y, and their derivatives at the point r: 

e'(r))^^^ = ao + aiX(r) + a2Y{v) + a,V^X{v) 



X,Y 



+ a,X\r) + a^{VXf{v) + ... , (51) 



:\v) + a,{yxf{ 

ri\r))^^ = bo + biX{r) + b2Y{r) + bsV^X{r) 



+ biX\r) + b^iyxyiv) + ... . (52) 

The terms on the right sides of (|51|) and (|52|) include all possible rotationally invariant 
combinations of X and Y and their derivatives. The constraint of rotational invariance 
arises from the rotational invariance of the free action (^). In (^Tj) and (0), we have 



shown explicitly only those terms that will ultimately be needed to calculate the quantum 
corrections to the density profile. 

The coefficients and bi in (^1]) and (^2]) can be reduced to integrals over a wave- vector 
k, as illustrated by the explicit calculation of a diagram presented in Appendix |B[ Having 
expanded the sources as Taylor series around r, the only scale in the integrand is A. By 
dimensional analysis, a convergent integral must have the form of the appropriate power 
of A multiplied by a numerical coefficient. However some of the integrals have infrared 
or ultraviolet divergences, and thus require infrared or ultraviolet cutoffs. The ultraviolet 
divergences either cancel in quantities such as p(r) and ('?/'(r)) or they are removed by 
renormalization. Infrared divergences reflect a failure of the assumption that the sources can 
be expanded in a Taylor series inside the loop integral. If these divergences do not cancel, it 
simply indicates a breakdown of the gradient expansion due to the sensitivity of the quantum 
corrections to nonlocal effects involving the length scale L. 

The propagator factor / dt'D^^yi^^ r', t') in (^) can be expanded in powers of X and its 
derivatives at the point r. The dependence on the source Y is removed by the integration 
over t', which corresponds to evaluating the Fourier-transformed propagator at a; = 0. Since 
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the off-diagonal components of tlie propagator P6|) vanisli at zero frequency, the source 
Y does not contribute. Examples of diagrams that contribute to / dt' D^y ^-^^ shown in 
Figure 8. The contribution from the first diagram is given by the upper- left component of 



the propagator matrix in (46) 



j dt' D^^{r,r',t') = -2m J 



d^k 
(2^ 



-jk-(r-r')_ 



1 



fc2 + A2 



(53) 



The other diagrams in Figure 8 involve integrals over the positions r" of sources X{r"). In 
coordinate space, the propagator factor ( p3D falls exponentially when |r — r'| exceeds 1/A. 
If we assume that the source X varies significantly only over a much greater length scale 
L, then we can expand X{r") as a Taylor series around the point r" = r. The function 
/ (it' D^y(r, r', t') can then be expressed in terms of X{r) and its derivatives at the point r. 
The terms coming from the diagrams in Figure 8 include 



r, r 



',t') 



-2m 



d^k 



(27r)3 
- 2iViX(r 



-jk-(r-r') 



+ X(t) 



1 



(P+A2)3 



P + A2 

+ V,V,X(r) 



(A;2 + A2)2 



(A;2 + A2)s 



+ 2V,XV,-X(r) 



6- 



(A;2 + A2)4 (A;2 + A2) 




The complete expression involves all possible powers of X and gradients of X. 

The result (^^ can be used to express {ip{r)) and p(r) as an expansion in powers of X and 
Y and their derivatives at the point (r). In the expression for (^(r)) in (pO|), the propagator 
factor is integrated against a function /(r') that depends on the sources T, Z, X, and Y. 
The integral can be evaluated by expanding /(r') as a Taylor series around the point r' = r. 
Using the expression (^^ for the propagator factor, we can evaluate the integral over r'. The 
resulting expression for the integral includes the terms 



dt'D%y{r,r',t'))f{r') 



-2mi 



^ ^X(r) + ^V^X(r) + 



A2 + A4' 



A6 



A8 



/(r) 



+ -VX(r) ■ V/(r) + ^VV(r) + ... 



(55) 



Applying this formula to the integral in (^Dj) and using the expansions (^Tj) and (^) for 
(^^) and {f]"^), we obtain an expansion for (^) in powers of X, Y, and their derivatives. 
Inserting the expansions for (^), (^2^^ a^^j ^^2^^ ^^^^^ jjg^ 



and (|38D , we obtain expansions for 
the condensate and the density in powers of X, Y, and their derivatives. 



4 One-loop calculation 

In this section, we calculate the one-loop quantum corrections to the density profile p(r) and 
to the condensate profile (^(r)) to second order in the gradient expansion. The appropriate 
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choice for the background field v is the mean field which satisfies (pS]): 



v{v) = 0o(r) . (56) 

The quantum fields ^ and 77 in ( p6D describe quantum fiuctuations around the mean field. 
The condensate profile (|37|) reduces to 

(^(r)) = 0o(r) + 4^(e(r)), (57) 



p{r) = 02(r) + y20o(r)(e(r)) + ^ (^'(r)) + l{v'{r)) - 6p. (58) 



while the number density (BSI) reduces to 



The fact that the mean field </)o satisfies the classical equation ( P3[ ) can be used to simplify the 
expressions (^OD-(^) for the sources. We can also drop the counterterms 5/i and 6g in the 
sources X, Y, and Z. These sources appear only in diagrams that are at least first order in 
the loop expansion. The counterterms appearing in these sources are therefore needed only 
to cancel ultraviolet divergences that arise at second order or higher in the loop expansion. 
Thus the sources can be simplified to 

r(r) = 5/i0o(r) - ^0^(r), (59) 
X(r) = - 2m#^(r) - ^(r) , (60) 

F(r) = -^(r), (61) 
Z(j) = -|0o(r). (62) 

The expressions (0) and (0) for the condensate and the density are nonlocal functional 
of the mean field 0o- If {'^i'"^)) p{^) ^i-re expanded in powers of the sources X and Y and 
their derivatives at the point r, the expansions include infinitely many terms. They can be 
reduced to local functionals of 0o by consistently truncating the expansions. We will reduce 
(^) and ( ^81) to local equations at a specific point tq by (a) choosing a specific value for 
the arbitrary parameter A and (b) truncating the equations at second order in the gradient 
expansion. Note that the source Y in ( |6lD is already second order in the gradient expansion. 
Thus if we truncate the equations at second order in the gradient expansion, we need only 
include terms up to first order in Y and we can omit all derivatives of Y. We also need only 
include terms up to first order in V^X and up to second order in VX. However, we still 
must include all possible powers of X. 

In order to reduce the expansions for (ip) and p to a finite number of terms, we will 
choose A so that X(ro) is second order in the gradient expansion at a specific point fq. If 
we evaluate (ip) and p at the point tq and then truncate them at second order in gradients 
of 00, the resulting expressions for {ip{ro)) and p(ro) are algebraic functions of (po, V0o, and 
V^0o evaluated at the point tq. Since we could have chosen any particular point for rg. 
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these algebraic relations must hold at any point r. The most convenient choice for A is the 
wavenumber that appears in the Bogoliubov dispersion relation (^2|) for a homogeneous gas 
with number density (pKrQ): 

= 2m(?02(ro). (63) 
The source X and its derivatives at the point tq then reduce to 



^(ro) 

VX(ro) 
V2X(ro) 



-Amg 0oV0o(ro) 



-Amg 



+ (V0o)' (ro). 



(64) 

(65) 
(66) 



We proceed to calculate the one-loop correction to the condensate {ip{r)), which is given 
by (^). Inserting ( [^ ) and (^) into (0) and using ( |53| ) to evaluate the integral over r', 
we obtain an expansion for (^) in powers of the sources and their derivatives. Inserting the 
expressions 

(e(ro)) 



), (pl|), ([H^), and (|B^-(|BBD for the sources at the point ro, we obtain 
-V2mg(j)Q\ 



3ao + ^0 _ Sfx Sg 2 
2 % + /° 



1 

A2 



6fi 
9 

9 , 



9ai + 3bi + 3a2 + &2 1 



1 

A4 



3ao + bo- 4^ ) — + 
3ao + 6o - 4 

-(3a3 + 63 - 12a4 - ^h) + 2(3a5 + 65)A 



A2 



+ (3a3 + 63 



(67) 



The coefficients and 6j are given in Appendix y. The coefficients and 60 are cubically 
ultraviolet divergent, while ai and 61 are linearly divergent. The divergences are cancelled 
by taking the counterterms 6fi and 6g to have the values 



6fi 
59 



127r2 
1 



47r2 



(m^iAuv) g . 



(68) 
(69) 



The counterterm 5g in ( p9| ) agrees with that obtained by expanding (|T?p to first order in 
mg A\jY. Using the results for and bi in Appendix |C| and the value of A given in ([63|), the 
condensate at the point Tq reduces to 



487r2 



{2mg) 



3/2, 



167r2 



+ 



2m(7 
/29 



'49 5 8m^02^ 
18 ^ 24 



^IR 



0g 



1 , 8mg(f), 
log 



16""° A2 



IR 



(70) 



where Air is an infrared cutoff. The logarithmic infrared divergences arise from the coeffi- 
cients 62, &3, and 65. The divergences indicate that the gradient expansion for the condensate 
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breaks down at second order. Thus we can obtain a local expression for the condensate only 
to leading order in the gradient expansion. Keeping only the first correction term in ([70|) , 
the result is 

(V^(r)) = 0o(r) [l - -^^{2mgf^,{v)\ . (71) 

We derived this equation at the point tq defined by our choice (|63D for the arbitrary parameter 
A. However, our final result for {iP{yq)) is an algebraic expression in terms of 0o(ro)- Since 
we could have chosen any particular point for fq, that algebraic expression must be valid at 
any point r. 

We next calculate the one- loop corrections to the density, which is given by ([5^). The 
expression for (^) at the point tq is given by (|67|). The corresponding expressions for (^^) 
and ijf) are obtained by inserting the expressions (^) and for the sources at the 

point ro into (|5TD and (|5^). The resulting expression for the density at the point Tq is 



p(ro) 



bp - 2 \ 

9 9 



+ 



3ao + &o - 4— 
9 , 



3ao + 6o - 4 



bp 



-2(a3 



9 
6a4 - 



— + (4ai + 6i + aa) + SogA^ 
A^ 

I ^ + (3ai + 6i) 
264)A2 + 4a5A^ 



(72) 



After using the expressions (|68|) and (|69|) for the counterterms bp and 5(7, the only remaining 
ultraviolet divergence is a cubic divergence that can be cancelled by choosing the density 
counterterm to be 

= I^A^^v (73) 

The infrared divergent coefficients 62? ^3? and 65 have cancelled in the expression (|72D for 
the number density. Thus the density has a well-defined gradient expansion through second 
order, in contrast to the condensate. Our final expression for the number density, including 
one-loop quantum corrections, is 



Mr) 




1 - 



67r^ 
1 

167r2 



(2m^)3/Vo(r) 



2'mg 



41V' 



9 



+ 



113 (Vc 



(74) 



We derived this equation at the specific point tq. However, our final expressions for p(ro) 
is an algebraic expression involving (pQ, V0o; and V^0o evaluated at the point vq. Since 
we could have chosen any particular point for Tq, these algebraic relations must hold at any 
point r. 

Combining ( |7TD and (|71D, we obtain a local expression for the condensate in terms of the 
density that is correct to leading order in the gradient expansion: 

1 - -l^{2mg)'/'^{T) . 



(^(r)) = yp(i 



487r2 



(75) 
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This agrees with a result obtained recently by Timmermans, Tommasini, and Huang in . 



The choice ( |63D for A is not unique. Any choice that makes X(ro) second order in the 
gradient expansion will be equally acceptable and must give the same final answer. For 
example, we could have chosen 

A2 = 2mgct>l{vo) + ^(ro) • (76) 

00 

In that case, (p^ ) would be replaced by X{tq) = 0. Following the effects of this change 
through the calculation, we find that the coefficient 4ai of V^0o/0o in (72) is replaced by 



3ai. However, the term — oq in ( |72[ ) depends on A, which is given in ([761). When this term is 
expanded in powers of gradients of 00; it generates additional terms proportional to V^0o/0o 
that precisely cancel the change in ([7^). Thus we recover the same final result (|7^). 

Note that the counterterms (|68D, (|6^), and ([73|) do not depend on the potential V{r). 
Thus the ultraviolet divergences in one-loop diagrams are removed by the same renormal- 
izations that are required for a homogeneous Bose gas. 



5 Self-consistent one-loop calculation 

In this Section, we present a self- consistent one-loop calculation of the equation for the density 
profile p(r) to second order in the gradient expansion. The calculation involves taking the 
equations for the density in the Hartree-Fock approximation and expanding them around the 
Thomas- Fermi limit. The result is the differential equation (H) that generalizes the Gross- 
Pitaevskii equation by taking into account the leading effects of quantum fluctuations. 

The self-consistent one-loop equations can be expressed as classical fleld equations for 
the one-loop effective action We describe briefly the diagrammatic representation of 

these equations. They correspond to summing all connected diagrams with arbitrarily many 
one-loop subdiagrams, but no subdiagrams with two or more loops. These diagrams have the 
structure of tree diagrams, with one-loop corrections added to the vertices and arbitrarily 
many one-loop corrections inserted into the propagators. These diagrams can be calculated 
using the perturbative framework developed in Section |^. The sum of all such diagrams is 
independent of the choice of the background fleld v{r) in (pB)) . However, the sum of all such 
diagrams can be greatly simplifled by choosing the background fleld v so that the ground- 
state expectation values of the quantum flelds ^ and r] vanish. This choice eliminates all 
one-particle-reducible diagrams which can be disconnected by cutting a single ^ or rj line. 
The only diagrams that remain are one-particle-irreducible diagrams. 

With the Cartesian parameterization (|5B|), the choice of the background fleld that sim- 
plifles self-consistent one-loop calculations is the condensate itself: 

viv) = (^(r)). (77) 

With this choice, the flelds ^ and 1] represent the quantum fluctuations around the ground- 
state expectation value of ip. Since v is real- valued, the expectation value of r] vanished 
automatically and the condition (^) can be written 

(e(r)) = 0. (78) 
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Thus the background field v must be chosen self-consistently so that the quantum fiuctuations 
around that background average to zero. We will refer to the equation ([78|) as the tadpole 
equation, because the one-loop quantum corrections to this equation correspond to Feynman 
diagrams like those in Figure 6 and 7 that look like tadpoles. Using the tadpole equation, 
the density (pSf ) reduces to 

p(r) = v'{r) + ^-{e{r)) + ^ (r/^(r)) - 5p. (79) 

The ground-state expectation values in ( |75D and ( [TUD are nonlocal functionals of the 
background v{r). Our strategy is to use the gradient expansion to reduce these functionals 
to local functions involving v{r) and its derivatives. The tadpole equation (|78D then reduces 
to an algebraic relation between v{r) and its derivatives, while (^) expresses p in terms of 
V and its derivatives. If we eliminate v from these two equations, we obtain an algebraic 
relation between p and its derivatives. This is the differential equation for p(r) that includes 
self-consistent corrections from one-loop quantum fiuctuations. 

To calculate the one-loop quantum corrections, we use the decomposition (^) of the 



action for quantum fiuctuations around a general background field v. The free part (^Sj) 
involves only the quantum fields ^ and t], but introduces an arbitrary scale A. The interaction 
part (|i7|) involves sources T, X, Y, and Z that are given in (|iO|)-(^3D. At one-loop order, 
the tadpole equation states that the expression (^) for (C,) vanishes, which implies 

= T(r) + lzir){e{r))^^^ + ^-Z{r) {v\r)) . (80) 
Similarly, the expression (|79D for the number density reduces at one-loop order to 

Pir) = v\r) + lieir))^^^ + - 6p. (81) 

The equations ( pO|) and (|8TD are integral equations whose solutions give the condensate and 
the density in the Hartree-Fock approximation. The quantities (C^)x,y and {rf)x,Y in ( P^D 
and (|81|) can be expanded in powers of X and Y and their derivatives using (^) and ([5^). 
Since these expansions include infinitely many terms, the equations ( pOD and ( |HT| ) can be 
reduced to local equations only by consistently truncating the expansions. We will reduce 
them to local equations at a specific point fq by (a) using the classical equations for p and 
V to simplify the expressions for the sources, (b) choosing a specific value for the arbitrary 
parameter A, and (c) truncating the equations at second order in the gradient expansion. 
We begin by simplifying the sources X, F, and Z in (^), (|4^), and (^) by using the 



classical equations T(r) = and p(r) = t>^(r). Since X, Y, and Z appear only in one-loop 
diagrams, any quantum corrections to the sources contribute only at second order in the 
quantum loop expansion. Using T = 0, we can eliminate the potential V from X and Y. 
Using V = ^/p, we can express X, Y, and Z in terms of p only. We can also simplify T by 
setting V = ^/p in the terms proportional to the counterterms 6p and 6g. Finally, we can 
drop the terms in X and Z that involve the counterterm Sg, since it is needed only to cancel 
ultraviolet divergences that arise at two loops or higher in the quantum loop expansion. 
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Thus, the expressions for the sources can be reduced to 

1 



T(r) 

X(r) 

F(r) 
Z(r) 



9 2( ^ 



f(r) + -^V^f(r) 
Zm 



/ — 

2mgp(r) — ^— (r) 



VP 



^ ^> 

-fyp(r). 



(82) 

(83) 

(84) 
(85) 



Note that the source Y in (0) is aheady second order in the gradient expansion. Thus if we 
truncate the equations at second order in the gradient expansion, we need only include terms 
up to first order in Y and we can omit all derivatives of Y. We also need only include terms 
up to first order in V^X and up to second order in VX. However, we still must include all 
possible powers of X. 

In order to reduce the expansions for ( |5DD and ( |STD to a finite number of terms, we choose 
A so that X(ro) is second order in the gradient expansion at a specific point Tq. A convenient 



choice for A is the wavenumber that appears in the Bogoliubov dispersion relation 
a homogeneous Bose gas with number density p(ro): 

A^ = 2mgp{ro) . 

With this choice for A, the source X and its derivatives at the point tq reduce to 



for 



(86) 



(87) 



VX(ro) 
V2x(ro) 



(89) 



VP 

-4mg [VP^^VP+i^VPy] M 

We proceed to determine the differential equation satisfied by p in the self-consistent 
one-loop approximation. This equation can be obtained by solving (|8TD for the condensate 
V in terms of the density p and its derivatives, and then eliminating v from the tadpole 
equation (|80D. If the tadpole equation is evaluated at the point fq and then truncated at 
second order in the gradient expansion, it reduces to 

3ao + 6o , Sai + bi ^. . 802 + 62 
+ ^ ^(ro) 



= T(ro) + Z{ro) 



Y{ro) 



+ ^^V^X(ro) 



+ ^^^(VX)2(ro 



The expression ([5^ ) for the source T involves v and V^f . Solving 
in the quantum corrections, we obtain 



f r 



v^(r) 



2v/p(i 



ao + 60 - 25p , ai + bi 

\ ^ — X(r) 



+ ^V^X(r) 



+ 



2 

04 + 64 



+ 



(90) 

for f (r) to first order 
02 -I- 62 



F(r) 



X^fr) 



+ 



as + 65 



fVX) 



+ 



(91) 
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We then substitute this expression for v into the source T(r) in (B2) and expand to first order 
in quantum fiuctuations. After calculating the derivative V^f appearing in T, we can set 
r = To and then truncate at second order in the gradient expansion. Using the expressions 
(15) and (|87D -(^9D for the sources and their derivatives, the expression for T(ro) reduces to 



T 



+ 



'ao + bo 



ao 



4 

60 - 25p 



25 p ^ ai + 61 



+ 



ai 




+ - Vp) 



2^ 2m ^ _ 



(92) 



where we have used the expression ( PUj) for A. The last term in ( P^ ) can be dropped, because 
it is proportional to the classical equation (||). Its effects are therefore of second order in the 
quantum loop expansion. This eliminates all occurrences of the potential V in the quantum 
corrections. Inserting the resulting expression for T(ro) into (pOD, the tadpole equation 
reduces to 



- - |p) 



+ 




^ ) 9Vp - yP^P 



A^ + asA' 



(93) 



Using the results for the coefficients and bi given in Appendix |C| and using (|86|) to set 
A^ = 2mgp{rQ), the equation for p reduces to 







(P 



1 



487r' 



:{2mg) 



3/2 



2 



Pv^(r) + ^^'v^(r 



17 
24m 



2^V'^{v) + {V^Y{v 



(94) 



We derived this equation at the point tq defined by our choice ( p6D for the arbitrary parameter 
A. However, our final result is an algebraic equation relating y/p and its derivatives at the 
point tq. Since we could have chosen any specific point for fq, this algebraic relation must 
hold at any point r. Using (^ to eliminate g in favor of a and using dimensional analysis to 
insert the appropriate factors of h into our equation (0), we obtain the differential equation 
(I) for the density profile. 

The term in ( p4|) can be obtained from previous work on the homogeneous Bose gas. 



(95) 



Differentiating the result (|3TD for the energy density with respect to ^/p, we obtain 



dS 



1 + 



{2mgf/^^ 
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Multiplying by — |, we reproduce the and terms in (p^. The ^^^^y^ (Vy/p)^ 
terms in ( p4|) are new results. 

As a check of the equation (0), we can verify that our one-loop expression for p given 
in (|7^ satisfies (|94D after expanding to first order in the quantum fluctuations. There is an 
important qualitative difference between the approximate solution ([7^ ) and the solution to 



the self-consistent equation (plf). The solution to (94) has the correct qualitative behavior 
even outside the condensate. In this region, the density is very small and only the terms in 
(0) that are linear in y/p are important. The equation therefore reduces to 

^ (/.-y(r))Vp(r) + ^V2y/p(r). (96) 



The quantum correction terms in (^4]) were calculated using a gradient expansion that is 
valid only inside the condensate. However, since these terms are all higher order in ^/p, 
their effects are negligible outside the condensate and it does no harm to include them. In 
contrast, the approximate solution (|7^ has the wrong qualitative behavior when 0o is small, 
because it is dominated by the V^0o and (V0o)^/0o terms. Thus that solution can only be 
used inside the condensate. 



6 Polar field parameterization 

The one-loop calculations in Sections ^ and ^ were carried out using the Cartesian param- 
eterization of the quantum field given in (0). There is nothing special about this parame- 
terization aside from its simplicity. Other field parameterizations should give the same final 
result for physical quantities. An example of an alternative field parameterization is the 
'polar parameterization: 

V^(r,t) = y^t;2(r) + a(r,t) exp(za(r, t)) . (97) 

The advantage of this parameterization is that it eliminates infrared divergences from indi- 
vidual Feynman diagrams that contribute to the number density. In this Section, we verify 
that the polar parameterization gives the same equation for the density profile. We also use 
this parameterization to show that the gradient expansion of the density breaks down at 
fourth order. 

With the polar parameterization (p7|) , the choice for the background field v that simplifies 
self-consistent one-loop quantum corrections is the one specified by the tadpole equation 

(a(r)) = 0. (98) 

The expression (^2p for the number density reduces to 

p(r) = v^{v) - bp. (99) 

Thus the choice of the background f (r) implied by the tadpole condition (^) is 

.;(r) = + + ... . (100) 
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Since this expression involves the ultraviolet divergent constant (5p, v has no simple physical 
interpretation. It is best regarded as a theoretical construct that should appear only in 
intermediate stages of a calculation. The simplicity of the expression (pQ]) for p comes at the 
expense of the expression for the condensate. Expanding (^) as a power series in a and a 
and taking the ground-state expectation value, we obtain 



(^(r)) = v{r) 



1 



(«'W) + .... 



(101) 



The expansion ( |101| ) includes infinitely many terms and we have written explicitly only those 
terms that contribute at one-loop order. The expectation values of operators involving four 
or more powers of cr or a contribute at two-loop order or higher. 

We begin our calculation by inserting the parameterization (^) into the action (pH) and 
expanding in powers of the quantum fields a and a: 



S[^] = S[v] + 



dt 
1 



d^r< -Ta + - iaa 



iVa) 



2 

9 + ^9 
4 



aa) — — (Va) 



+ 



4mf ^ J 



a 



2m ^ ' 



+ 



5mf ^ 



ra(Va)^ + 



{Vvf 



a' + 



,(102) 



where T is the external source given in (^Ol) . The parameterization (|97D leads to an infinite 
series of momentum-dependent interactions. We have dropped terms that are fourth and 
higher order in the quantum fields, since they do not contribute to the one-loop quantum 
corrections to the density profile. It is convenient to introduce an arbitrary parameter A into 
the action by rescaling the quantum fields as follows: 



a(r, t) 



A 



/mg 
/mg 



A 



e(r,t), 



(103) 
(104) 



After these rescalings, we separate the action into a free part and an interaction part as in 
(0). The free part is identical to ( ^51 ) and the interaction part is 



j dt j dv|V2 



\j2'mgv^ 



4m 4m 4m 



^/Imgi 



mv 



A 



(105) 



where X, Z, U, S, and W are external sources that depend on the background v. 
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Z{v) 
U{r) 
S{r) 
W{r) 



3/2 



3mf(r) \2mgv'^{r) 
A2 



2 ^ / 



V 



2m(yff 2(r) 
2mgv'^{v) 



A2 



3/2 



2mf(r) y2m5ff2(r 



(107) 
(108) 
(109) 
(110) 



The arbitrary parameter A, which was introduced through the rescahng of the fields, appears 
in both the free part of the action and the interactions. We will exploit the arbitrariness of 
this parameter to simplify the calculation of quantum corrections. 

The tadpole equation (^) can be written (^(r)) = 0. To first order in the quantum 
corrections, this equation implies that {A / y/2rngv^)T plus the sum of all one-loop tadpole 
diagrams vanishes. The one-loop diagrams include all possible insertions of the sources X, 
U, and S. The one-loop tadpole equation can be written succinctly in the form 
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\j2mgv- 
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1^ («^> 
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Arav A 
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111) 



where (. . .)x,u,s denotes the ground-state expectation value in the presence of the sources 
X, U, and S, but with no other interactions. The expectation values in (|111|) are nonlocal 
functionals of these sources. After Taylor-expanding the sources around the point r, these 
functional can be expanded in powers of X, U, and S and their derivatives at the point r: 



x,u,s 



x,u,s 



x,u,s 



Co + CiU + C2S + C3U'^ + C4US + C5S'^ + .. 

do + diX + d2 V^f/ + ds V^S 

+ di {yuf + d^vu ■VS + 4 ysf + ... 

eo + ei X + 62 V^U + 63 V^S 
+ 64 yuf + 65 Vf/ • V5 + 66 (V^)' + . . . 



(112) 

(113) 
(114) 



The terms on the right sides of ( |112|) - (|114| ) include all combinations of X, f/, and their 
derivatives that are allowed by rotational symmetry. We have written explicitly only those 
terms that are required to calculate the equation for the density through second order in the 
gradient expansion. 



The right side of the equation ( 111 ) is a nonlocal functional of v. When it is expanded 
in powers of X, f/, S*, and their derivatives, there are infinitely many terms. The equation 
can be reduced to a local one only by consistently truncating the expansions. We will 
reduce ( |111|) to a local equation at a specific point Yq by (a) choosing a specific value for 
the arbitrary parameter A, and (b) truncating the equations at second order in the gradient 
expansion. Since the sources X and Z in (p.06|) and (|107|) are already second order in the 
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gradient expansion, we need only keep terms of first order in X and Z and we can omit any 
derivatives of X. Moreover, we need only include terms that are first order in V^t/ and V^S* 
and terms up to second order in 'VU and VS*. However the equation ( |111| ) still includes 
all possible powers of U and S. We can reduce this equation to a finite number of terms at 
a specific point tq by choosing A so that the sources U and S vanish at the point fq. The 
required value is 

= 2mgv^{ro). (115) 
With this choice of A, the sources on the right sides of ( |112| )-( [TT^ ) reduce to 
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(117) 
(118) 

(119) 
(120) 
(121) 

(122) 



In the expression (|116|) for X(ro), we have dropped the term involving the counterterm bg^ 
since it is needed only to cancel ultraviolet divergences that arise at two loops or higher. 
With the choice ( |115| ) for A, the tadpole equation (|1 1 1|) simplifies at the point fq to 
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(123) 



After inserting the expansions ( |112D -( PT^ ) into the tadpole equation ( |123D and evaluating it 
at To, we can truncate it at second order in the gradient expansion. Using the expressions 
( |116|) - (|122|) for the sources, we obtain an algebraic equation relating t>, Vf , and V^f at the 
point fq. To express this equation in terms of p and its derivatives, we eliminate v using 
( llOOp . Since the term 5p/(2y/p) in ( |10CI|) is first order in quantum fluctuations, it is needed 
only in the term T in (|123[ ). The source T then becomes 



T 




Vp + 




2 / y/p 2m p 



■ (124) 



The last term in ( |124| ) can be dropped, because it is proportional to the classical equation 
(0). Its effects are therefore second order in the quantum loop expansion. After inserting 
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[2^) into ( |123|) , the tadpole equation reduces to 



" - " - f ") ^ + i^'^ 




(is + de 



(125) 



The coefficients q, di, and Cj are given in Appendix |C[ The coefficients are all infrared finite 
but ultraviolet divergent. The ultraviolet divergences from individual diagrams are more 
severe than those encountered with the Cartesian field parameterization used in Sections ^ 
and ^. The integrals do and Cq diverge as the fifth power of the ultraviolet cutoff, but they 
cancel in the combination d^ — cq. The remaining ultraviolet divergences are cancelled by 
the counterterms 6g, and 6p, whose values are given in (|68|), (|6^), and (|73|) . 

Using the expression for the coefficients given in the Appendix, the equation ( |125| ) reduces 

to 
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(126) 



After using ( |115| ) to set A^ = 2mgp{rQ), we reproduce the self-consistent one-loop equation 
(0) for the density profile. 

There has been a previous attempt to calculate the quantum corrections to the Gross- 
Pitaevskii equation [^. The authors used the polar field parameterization (|97D, with the 
background field f (r) equal to the mean field 0o(i")- They dropped all terms in the action 
that were third order and higher in a and a, and they also dropped second order terms that 
involved gradients of 0o or a. The only terms remaining in the action that contribute to the 
density profile are 



S[^P] = SIM + dt dh\-{ 



[aa — era) — 




(127) 



With such a drastic truncation of the action, the quantum corrections that they ultimately 
calculate are of no relevance to the problem of atoms in a trapping potential. This is evident 
from the fact that the Bogoliubov dispersion relation (|32D never enters into the quantum 
corrections that they calculate. Thus, their approach is incapable of reproducing the known 
results for a homogeneous Bose gas. 

A comparison of the calculation above with that presented in Section ^ demonstrates that 
the Cartesian field parametrization is more efficient than the polar field parameterization for 
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explicit calculations. With the polar field parameterization, one avoids infrared divergent 
integrals at intermediate stages of the calculation, but this advantage is compensated by the 
fact that the integrals are more severely ultraviolet divergent. The simplicity of the relation 
(p9|) between p and v is compensated by a tadpole equation that is more complicated 

than the corresponding equation (pOD in the Cartesian field parameterization. 

The advantage of the polar field parameterization is that it avoids cancellations of infrared 
divergences between different diagrams. This makes it easier to identify the sources of 
infrared divergences that are responsible for the breakdown of the gradient expansion. We 
will use this parameterization to show that the gradient expansion of the density breaks 
down at fourth order. The component of the propagator matrix (|46|) that is most infrared 
sensitive is D^"^ . For small loop momentum /c, the frequency uo in the loop scales like Ak/{2m) 
and D"^^ scales like 2m/fc^. The most infrared singular diagrams are those for which all the 
lines are rj lines. The term in the tadpole equation (|111|) that is most infrared sensitive is 
{{'Vr])^)x,u,s, because the operator (Vr])^ creates two r] fines. In the expansion ( |114| ) for 
that matrix element, the most infrared singular terms are those that involve only the source 
S, which couples to a pair of rj lines. The infrared behavior of the coefficient of a term in 
( |114|) that involves m factors of V and n factors of S can be determined by simple power 
counting. The integrand has a factor of l/k"^ for each of the n + 1 propagators. There is 
a factor of /c^ for the operator (Vr^)^ and a factor of k'^ for each insertion of S. Finally, 
dimensional analysis requires that each factor of V be compensated by a factor of 1/k in 
the integrand. Thus the integral must scale like 

diu [ £k k^ (^) ^ {k'^)'' (y) ~ fdkk^-"'. (128) 



,P/ ' ' \k. 

An infrared divergence can appear only if m > 4. Thus infrared divergences first appear in 
the tadpole equation when it is expanded to fourth order in the gradient expansion. 

Our explicit one-loop calculation of (ip) in ( |7DD showed that the gradient expansion for 
the condensate breaks down at second order. The analysis presented above shows that the 
gradient expansion for the density p does not break down until fourth order. The breakdown 
of the gradient expansion implies that the quantum corrections depend on nonlocal effects 
involving the length scale L for significant variations in p. While we have identified the 
orders at which the gradient expansions break down, we have not identified any deep reason 
for the gradient expansion of the density to be better behaved than that of the condensate. 



7 Implications for Present Traps 

In this Section, we estimate the magnitude of the effects of quantum fluctuations for Bose- 
Einstein condensates in existing magnetic traps. For numerical estimates, we will use pa- 
rameters characteristic of the sodium experiment 0, which has achieved the highest density 
condensates to date. In this experiment, N ^ 5 x 10^ sodium atoms were condensed in a 
trapping potential with a length scale i ^ 2 pm. The S-wave scattering length for sodium 
atoms is a 0.005 pm. The number density that was attained at the center of the condensate 
was p ~ 400 //im^. 

Baym and Pethick have presented a simple qualitative analysis of the solution to the 
Gross-Pitaevskii equation that allows one to determine how various quantities scale with the 
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number of atoms [Q. The qualitative behavior of the solution to the Gross-Pitaevskii 
equation depends crucially on a dimensionless parameter ( given by 

C ^ {^f . (129) 

where i is the length scale associated with significant variations in the potential V{r). For a 
harmonic oscillator potential, i = fi/ y/muj, where a; is the angular frequency of the harmonic 
oscillator. If ( is less than or of order 1, the size of the condensate is comparable to £ and 
the number density inside the trap scales like N/i^. If ( is much greater than 1, the size L 
of the condensate scales like L = (i and the density inside the trap scales like N/{(i)^. The 
condensates in existing magnetic traps are characterized by values of ( that are significantly 
greater than 1. (For the sodium experiment, ( ^ 13.) We will determine how the correction 
terms in our equation for the density profile scale with in the case C ^ 1- The expansion 
parameter -y/pa^ for the low-density expansion scales like (a/i. Although ( is large for 
existing traps, a/i is tiny and the product (a/i is small. (For the sodium experiment, 
(a/i ^ 0.03.) The modes that dominate the quantum corrections have wavelengths on the 
order of 1/^/pa, which scales like Since this is small compared to the length scale (i 
for significant variations in p(r), it is reasonable to expand the quantum corrections using 
the gradient expansion. The gradient expansion corresponds to an expansion in powers of 
the dimensionless quantity l/{yfpa(^i), which scales like Thus, inside the condensate, 

quantum corrections are suppressed hj (^a/i and corrections from second order in the gradient 
expansion are suppressed by l/C^. 

Outside the condensate, the density p rapidly approaches 0. The only terms in the 
equation (Q) that are important in this region are the terms that are linear in ^/p. The 
scale of the gradient is now set by the length scale i for significant variations in V^(r). In 
this region, the basic assumption underlying our calculation, that the quantum corrections 
are dominated by wavelengths of order 1/y/pa, breaks down completely. However, all the 
quantum corrections terms are higher order in y/p and therefore have a negligible effect on 
the solution outside the condensate. Thus it does no harm to include the quantum correction 
terms in (^ in the exterior region. 

The crossover region between the interior and exterior of the condensate can be character- 
ized by the fact that the V^y/p term and the p^/p terms become comparable in importance. 
The gradient expansion breaks down in this region. At the beginning of the crossover re- 
gion, V still scales like but the density has decreased to the point that p scales like 
a/{(i)^. Therefore the quantum loop expansion parameter y/pa^ scales like a/{(i). As long 
as this quantity is sufficiently small, the quantum correction terms in (0) are negligible. (In 
the sodium experiment, we have a/{(i) ~ 0.0002.) Thus it does no harm to include the 
quantum correction terms in (|]) in the cross-over region. We conclude that the differential 
equation (Q) can be used to calculate the density profile everywhere. 

We now give a quantitative estimate of the error from truncating the quantum loop ex- 
pansion at one-loop order. A simple estimate of the relative magnitude of the quantum cor- 



rections is the ratio of the correction term in (^ to the p^/p term, which is (32/3)y'pa^/7r. 
For the sodium experiment, this ratio is approximately 0.04 at the center of the condensate. 
This is small enough that quantum corrections can be treated as small perturbation to the 
mean-field approximation. Since y/pa^ scales like N^^^, the number of atoms in the trap 
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could be increased by many orders of magnitude and the condensate would still be within 
the perturbative region. 

A naive estimate of the relative magnitude of two-loop quantum corrections is the square 
of the magnitude of the one-loop quantum corrections. Their effects should therefore be 
negligible. One complication is that the two- loop correction depends not only on the S-wave 
scattering length a, but also on a second parameter that represents a pointlike contribution 
to the 3 — > 3 scattering amplitude |T^. If this parameter is anomalously large, the two- loop 
quantum corrections could be significantly larger than the naive estimate. 

We next quantify the errors from truncating the equation for the density at second order 
in the gradient expansion. A simple estimate of the relative magnitude of contributions 
from second order in the gradient expansion is the ratio of the V^y^ term to the term 
in the Gross-Pitaevskii equation (|^). Assuming that V^^p scales like a/p/(C^)^? the ratio 
is l/(8pa(C£)^). In the sodium experiment, this ratio is roughly 10~^ at the center of the 
condensate. These corrections are therefore negligible. The relative importance of quantum 
corrections that are second order in the gradient expansion increases as one approaches the 
edge of the condensate, where the gradient expansion breaks down. However they are still 
suppressed by a quantum loop factor of order aj {(,€). 

The gradient expansion for the density breaks down at fourth order, but the breakdown 
in only logarithmic in 8pa(C-^)^. We can estimate the magnitude of these corrections by 
taking the logarithms to be of order 1. These corrections are therefore suppressed by two 
powers of l/(8pa((^£)^). There is also an additional suppression factor of since terms 

of fourth order in the gradient expansion enter only through quantum corrections. Thus 
these corrections should be completely negligible. 

In this paper, we have developed a framework for calculating the dominant effects of 
quantum fluctuations in a Bose-Einstein condensate trapped by an external potential. Our 
method is based on a combination of the Hartree-Fock approximation and an expansion 
around the Thomas- Fermi limit. We have illustrated the method by calculating the self- 
consistent one-loop equation for the density profile to second order in the gradient expansion 
and the relation between the condensate and the density to zeroth order in the gradient 
expansion. It should be straightforward to use this method to calculate the effects of quantum 
fluctuations on other properties of the condensate at zero temperature, such as the spectrum 
of its collective excitations. It would also be interesting to extend the method to nonzero 
temperature, so that one could study how the effects of quantum fluctuations vary with 
temperature. 

This work was supported in part by the U. S. Department of Energy, Division of High 
Energy Physics, under Grant DE-FG02-91-ER40690. We thank T.-L. Ho for valuable dis- 
cussions. 

A Integrals 

In this Appendix, we give analytic expressions for the frequency integrals and the wavevector 
integrals that are required to calculate the one-loop quantum corrections to the density and 
the condensate. 

Since time-independent sources do not change the energy, the frequency integrals are 
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rather simple. They can be evaluated using contour integration. The specific integrals that 
are required are 



du^ 1 ^ , (-l)-l-3---(2n-3) 1 



ox - ^(-1)"+ ^ 7 -^r-^. (Al) 



du; ^ , (-l)-l-3---(2n-3) 1 



Time-independent sources that are inserted into a loop diagram change the wavevector 
k of the propagators in the loop. The gradient expansion corresponds to expanding the 
loop integral in powers of the wavevectors pj of the sources and of the external lines of the 
diagram. After averaging over integration angles, the k integrals that are required have the 
form 

I = i^^^fT^ .^3. 

''"''^ - J (27r)3 [kVWn^Y ' ^ ^ 

where m and n are integers. If m and n satisfy m — n < — |, this integral is ultraviolet 
convergent. If 2m — n> —3, the integral is infrared convergent. If it is both ultraviolet and 
infrared convergent, its value is 

W^ ^'—^y '"' a— ■ ™.5<„<2,„.3. ,A4) 

The ultraviolet-divergent integrals that are required are Im,n for n — m = — 1,0,1, which 
have power ultraviolet divergences. A convenient way to regularize the integral is to subtract 
pure powers of k from the integrand that will remove the ultraviolet divergence, and then to 
add those powers of k back in with an ultraviolet cutoff k < Auy. The regularized integral 
is then 

1 roc, / L2+2m-n m-n+1 / _n\ \ 

T = _/ — V I 2 ] A2ij,2+2m-2n-2i \ 

- 2772 70 '^H(P+A2)W2 [, )^ I" ) 

1 m-n+1 /n\ /-Ayy Q 

+ A E n^'' dkk'+'"^-'--'\ n<m+^. (A5) 



i=0 



The first integral in (^^) is convergent and is equal to the expression on the right side of 
). The only dependence on the ultraviolet cutoff comes from the remaining integrals in 
), and they give a polynomial in Auy. The ultraviolet divergent integrals that arise in 
our calculation are 



'n— l,n 



2 Aw + ^ tn^ ' A , (A6) 
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'"+i,n - ^ Auv - 3 AuvA + ^ AuvA + A j . {AS) 
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The infrared-divergent integrals that are required are Ijn,n for n — 2m + 3, which have 
logarithmic infrared divergences. The integrals can be regularized by imposing an infrared 
cutoff k > Air. In the limit Air <^ A, the regularized integral reduces to 
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'm,2m+3 



27r2 Jo 



dk 
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The specific integrals that are needed in our calculations are 
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Note that the integral /_2,-i is ultraviolet divergent as well as infrared divergent. 



B Explicit Calculation of a Diagram 



In this Appendix, we illustrate the calculation of one-loop diagrams that contribute to the 
quantum corrections to the density by calculating one diagram in detail. We consider the 
last diagram in Figure 3, which represents a contribution to the matrix element (C^)x,y 
involving two insertions of the source X. 

It is convenient to calculate the diagram in wavevector space, and then Fourier transform 
to get the diagram in coordinate space. The diagram involves an integral over the frequency 
uj and an integral over the wavevector k running around the loop. Letting the wavevectors 
of the two sources be pi and p2, the expression for the diagram is 
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(Bl) 



We have written the Feynman rules for each of the propagators and vertices in the loop in 
the order in which they appear as you go clockwise around the loop. There is a symmetry 
factor of I and the factor of 2 inside the integral is the Feynman rule for the operator 
There is an implied -l-iO"^ prescription in the denominator of each of the propagators. 

The first step in evaluating the diagram is to expand the integrand to second order in 
powers of the external wavevectors pi and p2. The expansion of the denominators has the 
form 

1 



K|k + P|) 
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We can now use the formula (^) to evaluate the integral over uj. This reduces the diagram 
to an integral over k. We can average over the angles of k by making the substitutions 
k'^y k^6^^ /3 and /c* — > 0. After simplifying the diagram, we find the terms 
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There are also terms proportional to {p\ + p2)^(Pi)^(P2) that we have dropped. They 
correspond to terms of the form XV^X, which first contribute to the density at fourth order 
in the gradient expansion. Expressing the integrals over k in ( |B3| ) in terms of the integrals 
Im,,n defined in Appendix 0, the expression reduces to 



-X(pi)X(p2)/3,5 

lo 
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(35/6,9 - 10/4,7 + 3/2,5) Pi ■ P2X(pi)X(p2) . 



(B4) 



After Fourier transforming to coordinate space, this becomes 



-h^X\v) + (35/6,9 -10/4,7 + 3/2,5) (VX)2(r). 



16 



192 



(B5) 



^From this expression, we can now read off the coefficients 04 and 05 in the expansion (^1| 
for {e)x.Y. 



C Coefficients 

In this Appendix, we express the coefficients that appear in the calculation of one-loop quan- 
tum corrections to the density in terms of the integrals Im n thai were defined in Appendix 

E 

We first list the coefficients that are used in Sections |^ and |^ to calculate the condensate 
and the density profile using the Cartesian field parameterization. The coefficients in the 
expansion (^Tj) for (^^)x,y are 

«o = lh,i, (CI) 

ai = (C2) 

02 = -^^,1) (C3) 

as = ^(-10/5,7 + 13/3,5-/1,3), (C4) 
3 

04 = T^4,5, (C5) 

lb 

as = (35/6,9 - 10/4,7 + 3/2,5) . (C6) 
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The coefficients in the expansion (|52D for (jf)x,Y are 

&o = (C7) 

hi = -\lo,i, (C8) 

b2 = \l-2,-i, (C9) 

h = ^(2/3,5- 3/1,3 , (CIO) 

h = (Cll) 

h = -Y^(/4,7 + 2/2,5 + 10,3) • (C12) 

We next hst the coefficients that are required in Section |] to calculate the density proffie 
using the polar field parameterization. The coefficients in the expansion (|112|) for {^'^)x,u,s 
are 

(C13) 
(CM) 
(C15) 
(C16) 
(C17) 
(C18) 

The coefficients in the expansion (|113|) for ((V^)^)x,c/,5 are 
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(C24) 
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The coefficients in the expansion (|114| ) for {{'Vri)'^)x,u,s are 
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Figure Caption 

Figure 1. Diagrammatic representation of the components of the 2x2 propagator matrix: 
(a) the diagonal propagator for ^, (b) the diagonal propagator for rj, (c) the off-diagonal 
propagator for ^ and rj. 

Figure 2. Diagrammatic representation of the interaction vertices associated with the sources 
T, X, Y, and Z and the four-point couplings. 

Figure 3. One- loop Feynman diagrams that contribute to (^^). 

Figure 4. One- loop Feynman diagrams that contribute to {rj'^) . 

Figure 5. Feynman diagrams contributing to {^) that involve the source T. 

Figure 6. One-loop Feynman diagrams contributing to (0 that involve a pair of ^ lines 
produced by the source Z. 

Figure 7. One-loop Feynman diagrams contributing to (^) that involve a pair of rj lines 
produced by the source Z. 

Figure 8. Feynman diagrams that contribute to the propagator -Dx,y zero frequency. 
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